cd "<INSERT DIRECTORY NAME>"

*Read in inflation data
insheet using "inflation.csv", clear
save inflation_data, replace

*Read in full income dataset
use cps_all2, clear
drop if hflag==1 /*Drop these because Medicare and Medicaid values not available for 3/8 sample in 2014 survey year*/

*Create state groupings consistent with 1967 groupings
rename statefip state
recode state (6 = 6 "CA") (9 = 9 "CT") (11 = 11 "DE") (12 = 12 "FL") (13 = 13 "GA") (17 = 17 "IL") (18 = 18 "IN") (21 = 21 "KY") (22 = 22 "LA") (24 = 24 "MD") (29 = 29 "MO") (34 = 34 "NJ") (36 = 36 "NY") (39 = 39 "OH") (41 = 41 "OR") (42 = 42 "PA") (47 = 47 "TN") (48 = 48 "TX") (54 = 54 "WV") (23 25 33 44 50 = 70 "ME-MA-NH-RI-VT") (26 55 = 71 "MI-WI") (27 19 = 72 "MN-IA") (31 38 46 20 = 73 "NE-ND-SD-KS") (10 51 = 74 "DE-VA") (37 45 = 75 "NC-SC") (1 28 = 76 "AL-MS") (5 40 = 77 "AR-OK") (4 35 8 = 78 "AZ-NM-CO") (16 56 49 30 32 = 79 "ID-WY-UT-MT-NV") (2 53 15 = 80 "AK-WA-HI"), gen(state_group)

*Create income measures
gen pre_tax_pre_trans = inctot
gen post_tax_pre_trans = pre_tax_pre_trans - fiitax - siitax - fica
gen post_tax_post_trans = post_tax_pre_trans + foodstamp + schllunch + houssub
gen post_tax_post_trans_hi = post_tax_post_trans + pmvcare + pmvcaid + emcontrb
gen post_tax_post_trans_hi_78 = post_tax_post_trans_hi
replace post_tax_post_trans_hi_78 = post_tax_post_trans_hi - siitax_1978 if year<=1978
gen post_tax_pre_trans_hi = post_tax_pre_trans + pmvcare + pmvcaid + emcontrb
gen pre_tax_post_trans_hi = post_tax_post_trans_hi + fiitax + siitax + fica

gen post_tax_post_trans_hi_WTP = post_tax_post_trans + emcontrb + 0.3 * (pmvcare + pmvcaid)
gen post_tax_post_trans_hi_WTPH = post_tax_pre_trans + foodstamp + schllunch + 0.7 * houssub + emcontrb + 0.3 * (pmvcare + pmvcaid)
gen post_tax_post_trans_hi_MS = min(post_tax_post_trans_hi, post_tax_post_trans * (3/2))
gen post_tax_post_trans_hi_MSA = min(post_tax_post_trans_hi, post_tax_post_trans * (6/5))

gen market_income = incwage + incbus + incfarm + incretir + incint + incsurv + incdivid + incrent + incchild + incalim + incasist + incrann + incpens + incidr + incaloth + incdrt + incother + emcontrb
gen welfare_income = incwelfr + incssi + foodstamp + schllunch + houssub + pmvcaid
gen socins_income = incunern + incss + incgov + incunemp + incwkcom + incvet + incdisab + inceduc + pmvcare
gen govt_income = welfare_income + socins_income
gen govt_income_post_tax = govt_income - fiitax - siitax - fica

local income_vars pre_tax_pre_trans post_tax_pre_trans post_tax_post_trans post_tax_post_trans_hi post_tax_pre_trans_hi pre_tax_post_trans_hi  post_tax_post_trans_hi_WTP post_tax_post_trans_hi_MS post_tax_post_trans_hi_78 market_income welfare_income socins_income govt_income govt_income_post_tax post_tax_post_trans_hi_WTPH post_tax_post_trans_hi_MSA

keep `income_vars' year serial povfam state state_group famsize hhsize personWeight relate age race  marst

foreach x in `income_vars' {
	bysort year serial: egen `x'_hh = sum(`x')
	bysort year serial povfam: egen `x'_fam = sum(`x')
	replace `x'_hh = `x'_hh/sqrt(hhsize)
	replace `x'_fam = `x'_fam/sqrt(famsize)
}


*Identify household types
gen working_age_ind = age>=18 & age<65 & relate==101 
gen married_ind = marst==1 & relate==101 
gen black_ind = race==200 & relate==101 
bysort year serial: egen working_age_head = max(working_age_ind)
bysort year serial: egen married_head = max(married_ind)
bysort year serial: egen black_head = max(black_ind)
gen working_age = age>=18 & age<65
gen married = marst==1
gen black = race==200

*Convert year from survey year to income year
replace year = year-1

merge m:1 year using inflation_data, nogen
drop if year<1963
*Drop 1969 due to data problems, interpolate this year
drop if year==1969
save initial_data3, replace

/*Program to get poverty trends anchored in 1963*/
capture program drop anchor_it
program define anchor_it
	/*Find base threshold*/
	use initial_data3, clear
	keep if year==`1'
	local p_star = 100*`2'/`3'
	_pctile `4' [aweight=personWeight], p(`p_star')
	scalar thres = `r(r1)' 
	use initial_data3, clear
	/*Get inflation in right terms*/
	gen `5'_`1'_temp = 0
	replace `5'_`1'_temp = `5' if year==`1'
	egen `5'_`1' = max(`5'_`1'_temp)
	replace `5' = `5'/`5'_`1'
	/*Get poverty rate in each year*/
	gen in_poverty = `4' < thres*`5'
	collapse in_poverty [aweight=personWeight], by(year)
	replace in_poverty = 100*in_poverty
	set obs `=_N+1'
	replace year=1969 if year==.
	sort year
	replace in_poverty = (in_poverty[_n-1] + in_poverty[_n+1])/2 if year==1969
	save `4'_`5', replace
end 

/*Program to get poverty trends anchored in given year, relative standard*/
capture program drop anchor_it_rel
program define anchor_it_rel
	/*Find base threshold*/
	use initial_data3, clear
	keep if year==`1'
	local p_star = 100*`2'/`3'
	_pctile `4' [aweight=personWeight], p(`p_star')
	scalar thres = `r(r1)' 
	_pctile `4' [aweight=personWeight], p(50)
	scalar med_income = `r(r1)'
	scalar pct_median = thres / med_income

	use initial_data3, clear
	sort year `4'
	bysort year: egen sum_weights = total(personWeight)
	bysort year: gen cum_weights = sum(personWeight)
	gen ordering = cum_weights / sum_weights
	gen median_income_temp=-9999
	replace median_income_temp = `4' if ordering[_n-1]<.5 & ordering>=.5
	bysort year: egen median_income = max(median_income_temp)
	gen thres = pct_median * median_income
	
	gen in_poverty = `4' < thres
	collapse in_poverty [aweight=personWeight], by(year)
	replace in_poverty = 100*in_poverty
	set obs `=_N+1'
	replace year=1969 if year==.
	sort year
	replace in_poverty = (in_poverty[_n-1] + in_poverty[_n+1])/2 if year==1969
	save rel_`1'_`4', replace
end 

/*Program to get poverty trends anchored in 1963 for race and age groups*/
capture program drop anchor_it_group
program define anchor_it_group
	/*Find base threshold*/
	use initial_data3, clear
	keep if year==`1'
	local p_star = 100*`2'/`3'
	_pctile `4' [aweight=personWeight], p(`p_star')
	scalar thres = `r(r1)' 
	use initial_data3, clear
	/*Get inflation in right terms*/
	gen `5'_`1'_temp = 0
	replace `5'_`1'_temp = `5' if year==`1'
	egen `5'_`1' = max(`5'_`1'_temp)
	replace `5' = `5'/`5'_`1'
	/*Get poverty rate in each year*/
	gen in_poverty = `4' < thres*`5'
	keep if black == `6' | black == `7'
	keep if age>=`8' & age<=`9'
	collapse in_poverty [aweight=personWeight], by(year)
	replace in_poverty = 100*in_poverty
	set obs `=_N+1'
	replace year=1969 if year==.
	sort year
	replace in_poverty = (in_poverty[_n-1] + in_poverty[_n+1])/2 if year==1969
	save `4'_`5', replace
end 

/**********************
****Poverty series*****
***********************/
*Absolute poverty series, anchored to 1963 (Figure 2, 3, 4, 5, 7)
anchor_it 1963 36436 187258 pre_tax_pre_trans_fam cpi_u
anchor_it 1963 36436 187258 pre_tax_pre_trans_hh cpi_u
anchor_it 1963 36436 187258 post_tax_pre_trans_hh cpi_u
anchor_it 1963 36436 187258 post_tax_post_trans_hh cpi_u
anchor_it 1963 36436 187258 post_tax_post_trans_hi_WTP_hh cpi_u
anchor_it 1963 36436 187258 post_tax_post_trans_hi_hh cpi_u
anchor_it 1963 36436 187258 post_tax_post_trans_hi_hh cpi_u_rs
anchor_it 1963 36436 187258 post_tax_post_trans_hi_hh pce
anchor_it 1963 36436 187258 post_tax_post_trans_hi_hh adj_cpi_u_rs

*Sensitivity Analysis (Table 1)
anchor_it 1963 36436 187258 post_tax_post_trans_hh pce
anchor_it 1963 36436 187258 post_tax_post_trans_hi_WTP_hh pce
anchor_it 1963 36436 187258 post_tax_post_trans_hi_MS_hh pce
anchor_it 1963 36436 187258 post_tax_pre_trans_hi_hh pce
anchor_it 1963 36436 187258 pre_tax_post_trans_hi_hh pce
anchor_it 1963 36436 187258 post_tax_post_trans_hi_fam pce

*Other sensitivity analyses
anchor_it 1963 36436 187258 post_tax_post_trans_hi_78_hh pce
anchor_it 1963 36436 187258 post_tax_post_trans_hi_WTPH_hh pce
anchor_it 1963 36436 187258 post_tax_post_trans_hi_MSA_hh pce

*Absolute poverty series, anchored to 1963, by demographic group (Figure 6)
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 0 0 0 99/*non-black*/
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 1 1 0 99/*black*/
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 0 1 0 99/*any race*/

anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 0 0 0 17/*non-black*/
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 1 1 0 17/*black*/
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 0 1 0 17/*any race*/

anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 0 0 18 64/*non-black*/
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 1 1 18 64/*black*/
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 0 1 18 64/*any race*/

anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 0 0 65 99/*non-black*/
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 1 1 65 99/*black*/
anchor_it_group 1963 36436 187258 post_tax_post_trans_hi_hh pce 0 1 65 99/*any race*/

*Relative poverty series, anchored to 1963 (Figure 2, 7)
anchor_it_rel 1963 36436 187258 post_tax_post_trans_hi_hh

*Relative poverty series, anchored to 2019 (Figure 11)
anchor_it_rel 2019 33984 324754 post_tax_post_trans_hi_hh

*Absolute poverty series, anchored to multiple years (Figure 11, Appendix Figure D1)
anchor_it 1973 22973 207621 post_tax_post_trans_hi_hh pce
anchor_it 1983 35303 231700 post_tax_post_trans_hi_hh pce
anchor_it 1993 39265 259278 post_tax_post_trans_hi_hh pce
anchor_it 2003 35861 287699 post_tax_post_trans_hi_hh pce
anchor_it 2013 46269 313096 post_tax_post_trans_hi_hh pce
anchor_it 2019 33984 324754 post_tax_post_trans_hi_hh pce

*Poverty in 1963 and 2019, by state (Appendix Figure C1)
use initial_data3, clear
keep if year==1963
local p_star = 100 * 36436 / 187258
_pctile post_tax_post_trans_hi_hh [aweight=personWeight], p(`p_star')
scalar thres = `r(r1)' 
use initial_data3, clear
/*Get inflation in right terms*/
gen pce_1963_temp = 0
replace pce_1963_temp = pce if year==1963
egen pce_1963 = max(pce_1963_temp)
replace pce = pce/pce_1963
/*Get poverty rate in each year*/
gen in_poverty = post_tax_post_trans_hi_hh < thres*pce
keep if year==1963 | year==2019
collapse in_poverty [aweight=personWeight], by(year state_group)
replace in_poverty = 100*in_poverty
reshape wide in_poverty, i(state_group) j(year)

/****************************
****Income Distribution******
*****************************/
*Income distribution 1963 vs. 2019 (Figure 10)
use initial_data3, clear
keep if year==1963 | year==2019
gen pce_2019_temp = 0
replace pce_2019_temp = pce if year==2019
egen pce_2019 = max(pce_2019_temp)
replace pce = pce_2019/pce
replace post_tax_post_trans_hi_hh = post_tax_post_trans_hi_hh * pce
/*median*/
bysort year: sum post_tax_post_trans_hi_hh [aweight=personWeight], detail/*Correct median*/
/*poverty threshold*/
local p_star = 100*36436/187258
_pctile post_tax_post_trans_hi_hh if year==1963 [aweight=personWeight], p(`p_star')
scalar thres = `r(r1)' 
_pctile post_tax_post_trans_hi_hh if year==1963 [aweight=personWeight], p(50)
scalar med_inc = `r(r1)'
scalar thes_med_ratio = thres / med_inc
_pctile post_tax_post_trans_hi_hh if year==2019 [aweight=personWeight], p(50)
scalar med_inc_2019 = `r(r1)'
scalar list/*Correct threshold*/
/*percent in poverty*/
gen in_poverty = post_tax_post_trans_hi_hh < thres
bysort year: sum in_poverty [aweight=personWeight]/*Get back correct poverty rates*/
replace post_tax_post_trans_hi_hh = ceil(post_tax_post_trans_hi_hh/1000)*1000
replace post_tax_post_trans_hi_hh = 0 if post_tax_post_trans_hi_hh<0
collapse (sum) personWeight, by(year post_tax_post_trans_hi_hh)
reshape wide personWeight, i(post_tax_post_trans_hi_hh) j(year)
egen sum_personWeight1963 = sum(personWeight1963)
egen sum_personWeight2019 = sum(personWeight2019)
replace personWeight1963 = personWeight1963 / sum_personWeight1963
replace personWeight2019 = personWeight2019 / sum_personWeight2019
twoway (line personWeight1963 post_tax_post_trans_hi_hh) (line personWeight2019 post_tax_post_trans_hi_hh) if post_tax_post_trans_hi_hh<=150000


/*****************
****Dependence****
******************/
*Programs to interpolate 1969 and multiply percent by 100
capture program drop add_1969 
program define add_1969
	set obs `=_N+1'
	replace year=1969 if year==.
	sort year
end

capture program drop interpololate 
program define interpololate	
	replace `1' = (`1'[_n-1] + `1'[_n+1])/2 if year==1969
	replace `1' = 100* `1'
end

*Percent with market income<50% by age and race (Figure 8)
use initial_data3, clear
gen age_group = ""
replace age_group = "00_17" if age>=0 & age<=17
replace age_group = "18_64" if age>=18 & age<=64
replace age_group = "65_99" if age>=65
gen dependent = market_income_hh < 0.5 * post_tax_post_trans_hi_hh
save tempy_1, replace

*By age group
use tempy_1, clear
collapse (mean) dependent [aweight=personWeight], by(year age_group)
reshape wide dependent, i(year) j(age_group) string
add_1969
interpololate dependent00_17 
interpololate dependent18_64 
interpololate dependent65_99

use tempy_1, clear
keep if black==0
collapse (mean) dependent [aweight=personWeight], by(year age_group)
reshape wide dependent, i(year) j(age_group) string
add_1969
interpololate dependent00_17 
interpololate dependent18_64 
interpololate dependent65_99

use tempy_1, clear
keep if black==1
collapse (mean) dependent [aweight=personWeight], by(year age_group)
reshape wide dependent, i(year) j(age_group) string
add_1969
interpololate dependent00_17 
interpololate dependent18_64 
interpololate dependent65_99

*Across all age groups
use tempy_1, clear
collapse (mean) dependent [aweight=personWeight], by(year black)
reshape wide dependent, i(year) j(black)
add_1969
interpololate dependent0
interpololate dependent1 

use tempy_1, clear
collapse (mean) dependent [aweight=personWeight], by(year)
add_1969
interpololate dependent


*Distribution of dependency (Figure 9)
capture program drop dep_dist
program define dep_dist

use initial_data3, clear
keep if age>=`1' & age<=`2'
keep if black==`3' | black==`4'
replace market_income_hh = 0 if market_income_hh<0
drop if post_tax_post_trans_hi_hh<0
gen market_share = market_income_hh / post_tax_post_trans_hi_hh
collapse (p10) market_share_10 = market_share (p25) market_share_25 = market_share (p50) market_share_50 = market_share (mean) market_share_mean = market_share [aweight=personWeight], by(year)
set obs `=_N+1'
replace year=1969 if year==.
sort year
foreach x in market_share_10 market_share_25 market_share_50 market_share_mean {
	replace `x' = (`x'[_n-1] + `x'[_n+1])/2 if year==1969
	replace `x' = 100 * `x'
}
end

dep_dist 0 99 0 1
dep_dist 0 99 1 1
dep_dist 0 99 0 0

dep_dist 18 64 0 1
dep_dist 18 64 1 1
dep_dist 18 64 0 0

*Dependence in 1967 and 2019, all ages, by state (Appendix Figure C2)
use initial_data3, clear
keep if year==1967 | year==2019
gen dependent = market_income_hh < 0.5 * post_tax_post_trans_hi_hh
collapse (mean) dependent [aweight=personWeight], by(year state_group)
reshape wide dependent, i(state_group) j(year)
replace dependent1967 = 100 * dependent1967
replace dependent2019 = 100 * dependent2019

*Dependence in 1967 and 2019, aged 18-64, by state (Appendix Figure C3)
use initial_data3, clear
keep if age>=18 & age<=64
keep if year==1967 | year==2019
gen dependent = market_income_hh < 0.5 * post_tax_post_trans_hi_hh
collapse (mean) dependent [aweight=personWeight], by(year state_group)
reshape wide dependent, i(state_group) j(year)
replace dependent1967 = 100 * dependent1967
replace dependent2019 = 100 * dependent2019


/**********************************
*****Dollars in survey vs admin****
**********************************/
*Appendix Table A1
use cps_all2, clear
drop if hflag==1 /*Drop these because Medicare and Medicaid values not available for 3/8 sample in 2014 survey year*/

foreach x in foodstamp schllunch houssub pmvcare pmvcaid emcontrb  {
	gen double `x'_fund = `x' * asecwt
}
collapse (sum) foodstamp_fund schllunch_fund houssub_fund pmvcare_fund pmvcaid_fund, by(year)

foreach x in foodstamp schllunch houssub pmvcare pmvcaid {
	format `x'_fund %20.0f
}
replace year = year - 1

/***********************************
*****Create SPM anchored to 1967****
***********************************/
*Create SPM series (Figure 7)
use pub_spm_all, clear
replace year = year - 1
gen threshold = SPMu_PovThreshold1_cen_Metadj / SPMu_Metadj_cen /*Removes geographic adjustment from threshold (https://static1.squarespace.com/static/610831a16c95260dbd68934a/t/62b1e2896ec23c5da63ad80c/1655825034054/SPM-public-use-data-documentation_06282021.pdf)*/
gen income = SPMu_Resources2 

*Scale thresholds such that poverty rate = 18.6% in 1967
gen income_pov = income / threshold
_pctile income_pov if year==1967 [aweight=asecwt], p(18.6)
scalar thres_mult = `r(r1)' 
replace threshold = threshold * thres_mult

*Calculate poverty rate series
gen in_poverty = income < threshold
collapse (mean) in_poverty [aweight=asecwt], by(year)
replace in_poverty = in_poverty * 100

